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An overview of advanced dynamical algorithms capable of spanning the widely disparate 
time scales that govern the decay of metastable phases in discrete spin models is pre- 
sented. The algorithms discussed include constrained transfer-matrix, Monte Carlo with 
Absorbing Markov Chains (MCAMC), and projective dynamics (PD) methods. The 
strengths and weaknesses of each of these algorithms are discussed, with particular em- 
phasis on identifying the parameter regimes (system size, temperature, and field) in 
which each algorithm works best. 

1. Introduction 

One of the most challenging hurdles to be overcome before the goal of computa- 
tional materials design can be realized is the problem of disparate time scales. The 
time scales range from microscopic (typically an inverse phonon frequency, roughly 
10~i2 s) to the time over which devices made from the materials must operate (years 
to centuries). The microscopic time scale is much smaller than the time between 
clock ticks on today's computers, and the times one desires to simulate are much 
too long to be achieved with traditional algorithms. Consequently, overcoming 
the hurdle of disparate timescales will require advanced dynamic algorithms which 
perform faster-than- real-time simulations. Here we briefly review three algorithms 
that allow such simulations for metastable decay in discrete spin models. Note that 
since we are interested in the dynamics of the problem, we cannot use advanced 
algorithms that change the dynamics (such as cluster algorithms, multicanonical al- 
gorithms, hybrid Monte Carlo algorithms, etc). Each advanced dynamic algorithm 
has its strengths and weaknesses, so which one should be used will depend not only 
on the individual problem to be simulated, but also on the parameter regime in 
which the problem is to be studied. 

We concentrate on applying these advanced dynamical algorithms to the case 
of metastable decay of the square-lattice Ising model with periodic boundary con- 



ditions. 



'his simple model has been used to study the metastable state of nanosca. 



e 



magnetsE and the thermal decay of information stored on magnetic recording medial. 
The Hamiltonian is 'H=^J^(^i j) <^i'^j^H^^cTi, with the Ising spins cr=±l. The 
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first sum is an interaction between nearest neighbors on a square LxL lattice, and 
the second sum is the interaction with the external field H. The temperature is 
taken to be below the critical temperature, rc«2.26- • - J. We start with all spins 
up (equal to +1) and H negative. In addition to the Hamiltonian, the dynamic 
must be specified: here the dynamic is randomly choosing one of the spins and 
deciding whether or not to flip it using a Metropolis flip probability^. One quantity 
we are interested in is the average lifetime, (r), before the system reaches some cut- 
off magnetization (taken to be zero) as it decays to the equilibrium state. Even in 
this simple model, which decays by homogeneous nuclcation and growth, there are 
different regimes where (t) has different functional formsB. There are four relevant 
length scales in the problem (far below Tc). Besides L and the lattice spacing, a, 
these are the critical droplet radius Rc and the typical distance Ro that a super- 
critical droplet (one with R>Rc) grows before it encounters another supercritical 
droplet. The critical droplet radius for circular droplets is Rc=(Joo{T)/2\H\msp 
where Uoo{T) is the surface tension along a primitive lattice vector when L^oo, 
and rrisp is the spontaneous magnetization. For the square lattice both Uoo (T) and 
msp are known exactly. The 'cross-over dynamical phase diagram' for metastable 
decay is shown in Fig. 1. 

For flelds that are too strong the droplet picture of decay is not valid, this is the 
Strong Field (SF) regimaJ. For weaker flelds the nucleation of multiple droplets (MD 
regime)u leads to decay of the metastable state via an Avrami mechanismd since 
Rc-^Ro-^L. A conservative estimate for the 'mean-fleld spinodal' that separates 
the SF and MD regimes is given by Rc=a/2, independent of L. For weak flelds and 
low temperatures a nucleation of a single droplet fSD regime) leads to decay of the 
metastable state, since Rc<^L<^Ro. Our estimateo of the 'dynamic spinodal' (where 
Ro^L) that gives the limit of the SD regime is where the standard deviation of the 
lifetime is equal to (t)/2. This cross-over depends on system size, and Fig. 1 shows 
simulation values using 10'^ escapes for L=10 and L=100. At |fl|>2 J and very low 
temperatures a single overturned spirtis the nucleating dropletia and the 'dynamic 
spinodal' for our dynamic is given bjfl (4— _ffDSp)/T'=|lii(i)— 0.82. Only the last 
term is an adjustable parameter, which has been chosen to give agreement for L 
between 8 and 240. For other dynamics, such as sequential updates, the functional 
form for i/DSp is differentQ. The 'coexistence regime'a, with L<Rc, is not shown in 
Fig. 1. 

As seen in Fig. 2, (r), with units of Monte Carlo steps per spin (MCSS), can 
vary over more than ten orders of magnitude. The remainder of the paper de- 
scribes briefly the algorithms that allow such simulations to be performed. The 
data points in Fig. 2 were taken with different techniques described in this paper, 
including the Monte Carlo with Absorbing Markov Chains (MCAMC) algorithm 
(with up to s—i) and the Projective Dynamics (PD) algorithm with and with- 
out a moving wall. With proper use of the algorithms, one obtains the correct 
lifetimes. Consequently, which points were taken with which algorithms is not in- 
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Fig. 1. The 'cross-over dynamical phase diagram' for metastable decay of the square-lattice 
Ising model. The dashed line corresponds to the 'mean-field spinodal'. The 'dynamic spinodal' is 
shown for L=10 and L=100. This is a heavy solid line for the analytical expression for \H\>2J, 
and dotted lines which connect the data points. Shown are the Strong Field (SF), Multi-Droplet 
(MD), and Single Droplet (SD) decay regimes. 
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J/T ■ J/T ■ 

Fig. 2. The average lifetimes, (r), in units of MCSS, for L=10 and _L=100 are shown as functions 

of inverse temperature for H=—0.75J . The same data are plotted in (a) and (b). The vertical 
dashed line is J /Tc- In (a) L'^(t) is plotted and the data points lie on top of each other in the 
SD regime. In (b) the same data points lie on top of each other when (r) is in the MD and SF 
regimes. The solid lines are from T— »0 predictions, as described in the text. 
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dicated in the figure. First note that the same data are plotted in parts (a) and 
(b) of Fig. 2. In the SD regime, Fig. 2(a), the lifetimes for different L lie on top of 
each other when L'^{t) is plotted. In the MD regime. Fig. 2(b), the values of (r) 
for different L lie on top of each other. The solid lines in Fig. 2 are T^O predic- 
tions for the exponential dependence valid for 2/3<|iJ|/J<l from low-temperature 
predicUonsli. In the SD regime! (r)=AsDexp[(24J-14|i7|)/T]/L2 and in the MD 
regimeO to (T)=AMDexp[(28J— 16|if|)/3T]. Only the prefactors are fit to the data. 
The SD prefactor Asd =0.186 is set to agree with the data at T—OAJ, while the 
MD prefactor Amd=0.563 is set to agree with the L=100 data point at T=J. 

2. Constrained Transfer Matrix 

In the 1960's Langeri used field-theoretical arguments to show that the nucleation 
rate, and hence the lifetime, can be related to the imaginary part of the analytically 
continued free energy density, Xm{!F). In particular, (T)oc[Im(J^)] . Numerically 
it is possible to obtain estimates for Xm(J-) by finding all eigenvectors and eigen- 
values of the 2^x2^ transfer matrix for an £xoo lattice. This is accomplished by 
introducing constrained entropy densitiesll2l in a manner analogous to the source en- 
tropy of a stationary ergodic Markov information source. The constrained transfer 
matrix (CTM) calculations obtain imaginary parts of the free energy density from 
the imaginary parts of the constrained entropy density. Although the transfer ma- 
trix does imt contain any explicit dynamics, the CTM method gives Xm{!F) values 
that agreell3 with droplet predictions for the functional form of (r). It also gives 
values for metastable magnetizations which agree with those obtained from Avrami 
analysis of MD growthB. Fig. 3 shows the stable and metastable magnetizations 
from CTM calculations at \H\=J/2 for various strip widths £. It illustrates some 
of the strengths and weaknesses of the CTM method. 

• CTM strengths 

o Can obtain Xm{J-) directly. 

o Can obtain metastable magnetization easily. 

• CTM weaknesses 

o Only feasible for small strip widths. This restricts applications to either 

mean-field like models or d—2 models with short-range interactions, 
o Must use very high precision in the computer calculations, 
o Interpretation subtle because of 'lobes' in transfer-matrix spectrum. 

3. Monte Carlo with Absorbing Markov Chains 

The first algorithm to obtain very large increases in simulated times for faithful 
dynamics of discrete spin models was put forward about 25 years ago by Bortz, 
Kalos, and LebowitzM. This n-fold way method is extremely efficient for low-T 
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Fig. 3. The absolute values of the stable (closed symbols) and metastable magnetizations (open 
symbols) at H=—0.5J for £=6,7,8,9. Note that the crossover between the SD and MD regimes 
is seen near T=1.75J. 
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simulations. The n-fold way algorithm has been rewrittenE3 for the discrete-time 
algorithm used in this paper, and this version is presented here. The underlying 

r 2 r 2 

dynamics corresponds to a 2 x2 Markov matrix. Since this matrix is too large 
to diagonalizdij, we can deal with it one piece at a time within our Monte Carlo 
simulation. If we look at an sxs submatrix T of the original Markov matrix, one 
has a matrix associated with the absorbing Markov chain (AMC) 

M = f J?^''' ^^''^ ) (3.1) 

\ -rtsxg J- sxs J 

where R takes into account the q ways to get out of the s states included in the 
transient matrix T. Here I is the identity matrix, and is the matrix with all 
zero elements. The initial vector (in this section we use the convention common 
in mathematics where operators act to the left) has zeroes everywhere with the 
exception of a one in the current spin configuration, which must be one of the s 
transient configurations. The random time m for exiting from the s transient states 
is the solution of the equation 

vjT^^,e<r<v^T-:~}e (3.2) 

where r is a uniformly distributed random number between and 1 , and e is the 
vector with all elements unity. Given that the s transient states are exited at time- 
step TO, the vector of (unnormalized) probabilities to exit to a particular one of the 
q states is given by 

^ (3.3) 

Once the s transient states have been exited, a new transient subspace is chosen 
and the algorithm is repeated. It is only Eqns. (3.2) and (3.3) which needJxi_be 
used to incorporate the methodology of AMC into a Monte Carlo simulationliJ't^. 

For the isotropic square-lattice Ising model in a field, each spin belongs to one 
of n=10 classes. The first 5 classes (1<?<5) have spin up and 5— z nearest-neighbor 
spins which are up. The last 5 classes (6<i<10) have spin down and 10— i nearest- 
neighbor spins which are up. Let Ci be the number of spins in class z, and let 
be the probability of flipping a spin in class i given that that particular spin was 
chosen for a spin- flip attempt. 

If s— 1, only the current spin configuration is in the transient subspace. This 
corresponds to the discrete-time 7i-fold way algorithm. In that case, for our square- 



lattice Ising model, the transient matrix is Tixi=l — L ^ Si=i CiPi^l — L ^Qi 
which defines Qw- Then Eq. (3.2) becomes to— 1< 



n(r)/ln(l-L-^(5io)<m. For 



TORiAi=— L^ln(r')/(5io which 



QxQ-^L^ the time of exiting can be approximated \y^. i 
leads to the continuous-time version of the rt-fold way algorithm. The probability 
of flipping one of the spins in class i is given from Eq. (3.3) by CiPi/QiQ, now 
independent of m. 

• MCAMC strengths 
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o Can give MANY orders of magnitude speed-up, so extremely large life- 
times can be simulated. 

o Extremely efficient at low temperatures and small L, i.e., in the SD 
regime. 

o Can have smaller errors in values for (r) compared with the same number 
of metastable escapes using normal simulations. This is because the 
exact average lifetime for exiting from the s transient states is given by 
iff {I — T)^^e, which can be incorporated into (r). 

o Can calculate higher moments of the escape time. 

o Provides a direct connection between discrete-time Monte Carlo simula- 
tions and the continuous time used in the original n-fold way algorithm. 

o Can he parallelized in a non-trivial fashion, which has been done for 
s=lE3. This allows one to efficiently use multi-processor scalable ma- 
chines to simulate discrete-event processes on large lattices for long times. 

• MCAMC weaknesses 

o Requires tabulation of Markov sub-matrices, which can be complicated 
for large L or small H. 

o Requires a small number of spins that evolve relatively quickly, so system 
stays in Markov sub-matrix for a large number of spin flip attempts. 

o Easiest to program efficiently for a small number of spin classes. 

o The parallelized s=l version looses some o£its efficiency at low temper- 
atures, compared with the serial algorithmll3. 



4. Projective Dynamics 

In order to show the similarities between the Projective Dynamics (PD) methodEUll 
and the Transition Matrix Monte Carlo (TMMC) methodN we use the TMMC 
notation in this section. The dynamics of the single-spin flip process acting on the 
2^ state vector P is given by the continuous-time Markov matrix F, 

^^ = J2r{a\aW,t). (4.1) 

{a'} 

The idea behind both the PD and TMMC methods is to coarse-grain this process 
by lumping together some of the 2 states according to some variable C. The 
probability of going from lumped state C to C" is given by 

^ ' o-e{cr|C} cr'e{cr|C'} 

where n{C') is the number of original states lumped into state C . This yields 

dP{C,t) 



dt {-^^ 

{C} 



W(C\C')P(C' ,t) (4.3) 



Advanced Dynamic Algorithms ■ ■ ■ 9 



where P(C, t) is the probabihty of being in lumped state C at time t. For the 
PD method we choose C=M , so we lump on the magnetization, and since we are 
interested in metastable decay we work at finite field and temperatures well below 
Tc- In the TMMCEa one lumps on the energy, so C—E, and one is interested in 
the critical behavior and works at H=0 and temperatures near IL. One interesting 
result is that the TMMC method yields dynamics much differenttj from single spin- 
flip dynamics, due to the allowed 'microcanonical' mixing within the lumped states. 
However, lumping on M does not allow spin flips to occur without jumping out of 
the current lumped state — and it has been proven that within statistics the correct 
value of (r) is obtainedE3. This means that memory effects due to the lumping in 
AI do not enter into the calculation of (t). For the PD case, the lumped matrix is 
tridiagonal and one has W{M\M')=s{M')6M+2,M'+g{M')dM-2,M' which defines 
the growing and shrinking rates of the stable phase as g{M)=J2^^^{ci) MPi and 
s(M)=J2i^e{ci) MPi- Here {ci)M is the average number of spins in class i during 
the simulated escape given that the configuration has magnetization M . 

The random walk starts at M=L^ and terminates at M=0. We define h{M) as 
the total time spent in M . If L is even, h{M)—0 for M odd, since these magneti- 
zations are not possible. Then 

It is also possible to introduce a moving wall in the magnetization to force 
the system out of the metastable state, and to measure g{n) and s{n) during 
this (hopefully quasi-static) process. The wall position is given by Mvja\\{t)—{L'^ + 
1)— Vwaiii- For a soft wall the normal Monte Carlo flip probability is multiplied by 
Pwaii^exp [— c (Mncw ^ -^^waii)] if thc wall at position Mwaii is past the magnetiza- 
tion of the Monte Carlo new trial configuration. Another possibility is to introduce 
a hard wall, c=oo, that always flips a chosen up spin if Mwaii is past the current 
magnetization. Fig. 4 shows an example of lifetimes for different forcing walls as 
functions of the wall speed. The difference between the lifetimes measured from 
Eq. (4.4) and the number of MCSS before escape with the wall gives an estimate 
for the speed-up due to incorporating the forcing wall. Of course, for a wall moving 
too fast the process is not quasi-static and thc lifetimes are not accurate. How- 
ever, good results are obtained even for relatively fast walls. For example for the 
hard wall a velocity of 3xlO~'*M/L^MCSS gives a lifetime comparable to an actual 
escape but requires almost two orders of magnitude less computer time. The soft 
walls do not seem to help the convergence, and additional calculations are needed to 
calculate Pwaii ■ This makes the hard wall more computationally efficient compared 
with the soft walls. 



PD strengths 

o Can give MANY orders of magnitude speed-up, so extremely large life- 
times are possible. 
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Fig. 4. The average lifetime, (r), in units of MCSS is shown as a function of the velocity of the 
forcing wall (in units of M/L^MCSS). This is for a 10x10 lattice with T=0.9J and H=-0.75J. 
The data are for different degrees of softness of the wall, c=0.01 (□), c=1.0 (O), and a hard wall 
c=oo (v)- The filled symbols are (t) values from the growing and shrinking probabilities, while 
the corresponding open symbols are the measured number of MCSS to escape from the metastable 
state with the forcing wall. The dotted lines connect the data points and serve as guides to the 
eye. The dashed horizontal line corresponds to the result for (t) from standard simulations. 
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o Easy to obtain magnetization of metastable stale 
o Easy to obtain magnetization of saddle poinlO'El. 
o Allows extrapolation to very long lifetimes for weak IrEi. 
o Allows extrapolation to large L from simulations at smaller £E0. 
o Gives EXACT value of lifetime, except for statistical deviationst2l. 
o Has decreased errors in values for (t) compared with same number of 

metastable escapes using normal simulations, 
o Works well in both SD and MD regimes, as well as the coexistence regini£. 
o Can incorporate a 'moving forcing wall' to further enhance performanceES. 

PD weaknesses 

o Gives only approximate higher moments of lifetime due to memory effects 

ignored when the states are lumpedEJ. 
o The origiiial Markov chain is not 'lumpable', nor is it even 'weakly 

lumpableila. 

o It is nontrivial to know how fast to move the wall. 



5. Discussion and Conclusions 

We have demonstrated several dynamically faithful algorithms to study homoge- 
neous nucleation and growth in discrete spin models. These include the Constrained 
Transfer Matrix (CTM) method, the Monte Carlo with Absorbing Markov Chains 
(MCAMC) method which is a generalization of the n-fold way algorithm, and the 
Projective Dynamics (PD) method. Each method has its own strengths and weak- 
nesses, which are detailed in Sec. 2, 3, and 4, respectively. Consequently, which 
method should be used depends on the model to be studied and on the param- 
eter regime where the model is to be studied. Application and generalization of 
these methods to faithfully study the dynamics of other systems, including_SYSj- 
tems with randomness, continuous spin modelsEa, higher-dimensional systemalU'ta, 
and systems evolving via stochastic differential equations are some possible future 
interesting applications. 
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